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ABSTRACT 


Very  large  »eale  matrix  problems  currently  arise  in  the  context 
of  accurately  computing  the  coordinates  of  points  on  the  surface  of 
the  earth.  Here  geodesists  adjust  the  approximate  values  of  these 
coordinates  by  computing  least  squares  solutions  to  large  sparse  systems 
of  equations  which  result  from  relating  the  coordinates  to  certain  ob¬ 
servations  such  as  distances  or  angles  between  points.  The  purpose  of 
this  paper  is  to  suggest  an  alternative  to  the  formation  and  solution  of 
the  normal  equations  for  these  least  squares  adjustment  problems.  In 
particular,  it  is  shown  how  a  block -orthogonal  decomposition  method  can 
be  used  in  conjunction  with  a  nested  dissection  scheme  to  produce  an 
algorithm  for  solving  such  problems  which  combines  efficient  data 
management  with  numerical  stability.  As  an  indication  of  the  magnitude 
that  these  least  squares  adjustment  problems  can  sometimes  attain,  the 
forthcoming  readjustment  of  the  North  American  Datum  in  1985  by  the 
National  Geodetic  Survey  is  discussed.  Here  it  becomes  necessary  to 
linearize  and  solve  an  overdetermined  system  of  approximately  6,000,000 
equations  in  V00,000  unknowns  -  a  truly  large-scale  matrix  problem. 
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Introduction. 


Recent  technological  advances  have  made  possible  the  collection  of 
vast  amounts  of  raw  data  describing  certain  physical  phenomena.  As  a 
result ,  the  sheer  volume  of  the  data  has  necessitated  the  development 
of  new  elaborate  schemes  for  processing  and  interpreting  it  in  detail. 

An  example  is  in  the  adjustment  of  geodetic  data. 

Geodesy  is  the  branch  of  applied  mathematics  which  is  concerned 
with  the  determination  of  the  size  and  shape  of  the  earth,  the  directions 
of  lines  and  the  coordinates  of  stations  or  points  on  the  earth's  surface. 
Applications  of  this  science  include  mapping  and  charting,  missile  and 
space  operations ,  earthquake  prediction, and  navigation.  The  current  use 
of  electronic  distance  measuring  equipment  and  one-second  theodolites 
for  angle  measurements  by  almost  all  surveyors  necessitates  modern  ad¬ 
justment  procedures  to  guard  against  the  possibility  of  blunders  as  well  as 
to  obtain  a  better  estimate  of  the  unknown  quantities  being  measured.  The 
number  of  observations  is  always  larger  than  the  minimum  required  to 
determine  the  unknowns.  The  relationships  among  the  unknown  quantities 
and  the  observations  lead  to  an  overdetermined  system  of  nonlinear  equations 
The  measurements  are  then  usually  adjusted  in  the  sense  of  least  squares 
by  computing  the  least  squares  solution  to  a  linearized  form  of  the  system 
that  is  not  rank  deficient. 

In  general,  a  geodetical  position  network  is  a  mathematical  model 
consisting  of  several  mesh-points  or  geodetic  stations,  with  unknown  posi¬ 
tions  over  a  reference  surface  or  in  three-dimensional  space.  These  stations 
are  normally  connected  by  lines,  each  representing  one  or  more  observations 
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involving  the  two  stations  terminating  the  line.  The  observations  may  be 
angles  or  distances, and  thus  they  lead  to  nonlinear  equations  involving, 
for  example,  trigonometric  identities  and  distance  formulas  relating 
the  unknown  coordinates.  Each  equation  typically  involves  only  a  small 
number  of  unknowns. 

As  an  illustration  of  the  sheer  magnitude  that  some  of  these  problems 
can  attain,  we  mention  the  readjustment  of  the  North  American  Datum  - 
a  network  of  reference  points  on  the  North  American  continent  whose 
longitudes,  latitudes  and,  in  some  cases,  altitudes  must  be  known  to  an 
accuracy  of  a  few  centimeters.  This  ten-year  project  by  the  U.S.  National 
Geodetic  Survey  is  expected  to  be  completed  by  198j.  The  readjusted  net¬ 
work  with  very  accurate  coordinates  is  necessary  to  regional  planners, 
engineers  and  surveyors,  who  need  accurate  refer  nee  points  to  make  maps 
and  specify  boundary  lines;  to  navigators;  to  oad  builders;  and  to  energy 
resource  developers  and  distributors.  Very  briefly,  the  problem  is  to  use 
some  6,000,000  observations  relating  the  positions  of  approximately 
200,000  stations  (400,000  unknowns)  in  order  to  readjust  the  tabulated  values 
for  their  latitudes  and  longitudes.  This  leads  to  one  of  the  largest  single 
computational  efforts  ever  attempted  -  that  of  computing  a  least  squares 
solution  of  a  very  sparse  system  of  6,000,000  nonlinear  equations  in 
400,000  unknowns.  This  problem  is  described  in  detail  by  Meissl  [I979l> 
by  Avila  and  Tomlin  [19791>  and  from  a  layman's  point  of  view  by  Kolata 
[1978]  in  Science. 

In  general  then,  geodetical  network  adjustment  problems  can  lead 
(after  linearization)  to  a  very  large  sparse  overdetermined  system  of  m 
linear  equations  in  n  unknowns 


t  f 


Ax  s*  b 


(l.l) 


where  the  matrix  A  ,  called  the  observation  matrix,  has  full  column 
rank.  The  least  squares  solution  to  (l.l)  is  then  the  unique  solution 


to  the  problem: 


m£n||b-  Ax||s 


An  equivalent  formulation  of  the  problem  is  the  following:  one  seeks  to 
determine  vectors  y  and  r  such  that  r  +  Ay  =  b  and  A^r  =  0 
The  least  squares  solution  to  (l.l)  is  then  the  unique  solution  y  to 
the  nonsingular  system  of  normal  equations 


A^Ay  =  aH> 


(1.2) 


The  linear  system  of  equations  (1.2)  is  usually  solved  by  computing 


the  Cholesky  factorization 

AtA  ~  R^R  ,  R  =  ^ 

• 

and  then  solving  Rw  =  A^b  by  forward  substitution  and  Ry  =  w  by 
back  substitution.  The  upper  triangular  matrix  R  is  called  the 
Cholesky  factor  of  A  . 

Most  algorithms  for  solving  geodetic  least  squares  adjustment  problems 
(see  Ashkenazi  [1971],  Bomford  [1971],  Meissl  [1979]  or  Avila  and  Tomlin 
[1979]) typically  involve  the  formation  and  solution  of  some  (weighted) 
form  of  the  normal  equations  (l.2)>  But  because  of  the  size  of  these 
problems  and  the  high  degree  of  accuracy  desired  in  the  coordinates,  it 
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is  important  that  particular  attention  be  paid  to  sparsity  considerations 
when  forming  A^A  as  well  as  to  the  numerical  stability  of  the  algorithm 
being  used.  It  is  generally  agreed  in  modern  numerical  analysis  theory 
(see  Golub  [1965],  Lawson  and  Hanson  [1974]  or  Stewart  [ 1978]  )that  ortho¬ 
gonal  decomposition  methods  applied  directly  to  the  matrix  A  in  (l.l)  are 
preferable  to  the  calculation  of  the  normal  equations  whenever  numerical 
stability  is  important.  Since  A  has  full  column  rank,  the  Cholesky 
factor,  R  ,  of  A  can  be  computed  by 


A  = 


L°J 


,  Q  Q  =  I  ,  R  = 


(1.3) 


where  the  orthogonal  matrix  Q  results  from  a  finite  sequence  of 
orthogonal  transformations,  such  as  Householder  reflections  or  Givens 
rotations,  chosen  to  reduce  A  to  upper  triangular  form. 


Since  A  has  the  orthogonal  decomposition 


Q 


then  defining  Q,  b  = 


,  where  c  is  an  n  -  vector. 


the  least  squares  solution  y  to  (l.l)  is  obtained  by  solving  Ry  =  c 
by  back  substitution.  The  greater  numerical  stability  of  the  orthogonal 
decomposition  method  results  from  the  fact  that  the  spectral  condition 
number  of  A^A  in  the  normal  equations  (l. 2)  is  the  square  of  the  spectral 
condition  number  of  A  .  The. orthogonal  decomposition  method  (l.j)  has 
other  advantages,  including  the  ease  with  which  updating  and  downdating  of 
the  system  (l.l)  can  be  accomplished,  and  the  fact  that  possible  fill-in 
in  forming  the  normal  equations  is  avoided  (see,  for  example,  BjBrck  [1976]) 
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However ,  orthogonal  decomposition  techniques  for  solving  large  sparse 
least  squares  problems  such  as  those  in  geodesy  have  generally  been 
avoided,  in  part  because  of  tradition  and  in  part  because  of  the  lack 
of  effective  means  for  preserving  sparsity  and  for  managing  the 
data. 

Modern  techniques  for  solving  large  scale  geodetic  adjustment 
problems  have  involved  the  use  of  a  natural  form  of  nested  dissection, 
called  Helmert  blocking  by  geodesists,  to  partition  and  solve  the  normal 
equations  (l.2).  Such  techniques  are  described  in  detail  in  Avila  and 
Tomlin  [l979]>  in  Hanson  [1978],  and  in  Meissl  [1979]  where  error  analyses 
are  given. 

The  purpose  of  this  paper  is  to  develop  an  alternative  to  the  formation 
and  solution  of  the  normal  equations  in  geodetic  adjustments.  We  show  how 
the  orthogonal  decomposition  method  can  be  combined  with  a  nested  dissection 
scheme  to  produce  an  algorithm  for  solving  such  problems  that  combines 
efficient  data  management  with  numerical  stability. 

In  subsequent  sections  the  adjustment  problem  is  formulated,  and  it 
is  shown  how  nested  dissection  leads  to  an  observation  matrix  A  in  (l.l) 
of  the  special  partitioned  form 


where  the  diagonal  blocks  are  normally  rectangular  and  dense  and  where 
the  large  block  on  the  right-hand  side  is  normally  sparse  with  a  very 
special  structure.  The  form  (l.k)  is  analyzed  and  a  block- orthogonal, 
decomposition  scheme  is  described.  The  final  section  contains  some 
remarks  on  the  advantages  of  the  approach  given  in  this  paper  and 
relates  the  concepts  mentioned  here  to  further  applications.  Numerical 
experiments  and  comparisons  are  given  elsewhere  in  Golub  and  Plemmons 
[1980]. 

2.  Geodetic  Adjustments. 

In  this  paper  we  consider  geodetical  position  networks  consisting 
of  mesh-points,  called  stations,  on  a  two-dimensional  reference  surface. 
Associated  with  each  station  there  are  two  coordinates.  A  line  connecting 
two  stations  is  roughly  used  to  indicate  that  the  coordinates  are  coupled 
by  one  or  more  physical  observations.  Thus  the  coordinates  are  related 
in  some  equation  that  may  involve,  for  example,  distance  formulas  or 
trigonometric  identities  relating  angle  observations.  An  example  of  such 
a  network  appears  in  Figure  1. 


FIGURE  1 


A  15  station  network. 


* 


mmmmm 


More  precisely,  one  considers  a  coordinate  system  for  the  earth 
and  seeks  to  locate  the  stations  exactly,  relative  to  that  system. 

Usually  coordinates  are  chosen  from  a  rectangular  geocentric  system  (see 
Bomford  [19713).  Furthermore,  a  reference  ellipsoid  of  revolution  is 
chosen  in  this  set  of  coordinates  and  the  projection  of  each  station  onto 
this  ellipsoid  determines  the  latitude  and  longitude  of  that  station. 

As  indicated  initially  in  Section  1,  the  relationships  among  the 
coordinates  of  the  stations  in  the  geodetic  network  lead  to  an  over¬ 
determined  system  of  nonlinear  equations 

'  F (p)  =  q  (2.1) 


where 

p  =  vector  of  unknown  coordinates,  and 
q  =  vector  of  observations. 

The  components  of  F(p)  represent  the  equations  that  express  the  relation¬ 
ships  among  the  unknown  parameters  and  the  observations  or  measurements 
made,  for  example,  by  surveyors. 

A  common  procedure  for  solving  the  overdetermined  system  (2.l)  is  the 
method  of  variation  of  parameters.  (This  is  generally  called  the  Gauss- 
Newton  nonlinear  least  squares  algorithm  in  the  mathematical  literature). 
Approximate  coordinates  are  known  a  priori.  Let 

p®  =  current  vector  of  approximate  coordinates. 

Then  if  F  has  a  Taylor's  series  expansion  about  p°  ,  there  follows  the 
relationship 
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F(p)  =  F(p°)  +  F'(p°)(p  -  p°)  +  ... 

where  F* (p^)  denotes  the  Jacobian  of  F  at  p°  .  Then  taking 

A  =  F* (p°) 

0 

x  =  p  -  p 
b  =  q  -  F(p^) 

and  truncating  the  series  after  2  terms ,  one  seeks  tbe  solution  to: 

min||b  -  Ax|L  .  (2.2) 

x  * 

The  least  squares  solution  y  then  represents  the  correction  to 
p°  .  That  is,  one  takes 

1 

P  =  P  +  y 

as  the  next  approximation  to  p  .  The  process  is,  of  course,  iterative 
and  one  can  use  p1  to  compute  a  further  approximation  to  p  .  Normally, 
the  initial  coordinates  have  sufficient  accuracy  for  convergence  of  the 
method,  but  the  number  of  iterations  is  often  limited  by  the  sheer  magnitude 
of  the  computations.  Thus  a  very  accurate  approximation  to  y  is  desired. 

Actually,  the  equations  are  usually  weighted  by  use  of  some  positive 
diagonal  matrix  W  ,  where  the  weights  are  chosen  to  reflect  the  confidence 
in  the  observations:  thus  (2.2)  becomes 

min||W^b  -  W^AxL  . 
x  *= 

For  simplicity,  we  will  use  (2.2)  in  the  analysis  to  follow.  The  procedure 
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we  discuss,  however,  will  not  be  complicated  by  the  weights. 

Due  to  the  sheer  volume  of  the  data  to  be  processed  in  many 
adjustment  problems,  it  is  imperative  to  organize  the  data  in  such  a 
way  that  the  problem  can  be  broken  down  into  meaningful  mathematical 
subproblems  which  are  connected  in  a  well-defined  way.  The  total 
problem  is  then  attacked  by  "solving"  the  subproblems  in  a  topological 
sequence.  This  "sub structuring"  or  "dissection"  process  has  been 
used  by  geodesists  for  almost  a  century.  The  method  they  have  employed 
dates  back  to  Helmert  [l88o]  and  is  known  as  Helmert  blocking  (see 
Wolf  [1978]  for  a  historical  discussion). 

In  Helmert  blocking,  geographical  boundaries  for  the  region  in 
question  are  chosen  to  partition  it  into  regional  blocks.  This  technique 
orders  the  stations  appropriately  in  order  to  establish  barriers  which 
divide  the  network  into  blocks.  The  barriers  are  chosen  so  that  the 
interior  stations  in  one  block  are  not  coupled  by  observations  to  interior 
stations  in  any  other  block.  These  interior  blocks  are  separated  by  sets 
of  junction  stations  which  are  coupled  by  observations  to  stations  in  more 
than  one  block.  An  example  of  such  a  partitioning  of  the  geodetic  network 
in  Figure  1  to  one  level  of  Helmert  blocking  is  provided  in  Figure  2. 

Here  the  circled  nodes  represent  the  junction  stations  chosen  for  this 
example . 

I 

/  *  • 
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1  '  FIGURE  2 

,  One  level  of  Helmert  blocking. 

»  The  particular  form  of  Helmert  blocking  we  will  use  here  is  the  same 

as  that  used  by  Avila  and  Tomlin  [1979]  for  partitioning  the  normal 
equations.  That  procedure,  in  certain  respects,  is  a  variation  of  the 
nested  dissection  method  developed  by  George  [197?],  [1977]; 

«• 

George  and  Lui  [1978] ;and  George,  Poole  and  Voight  [1978].  The  primary 
emphasis  of  the  nested  dissection  strategy  has  been  on  solving  symmetric 
positive- definite  systems  of  linear  equations  associated  with  finite  element 
!  schemes  for  partial  differential  equations.  There,  the  finite  element  nodes 

are  ordered  in  such  a  way  that  the  element  matrix  B  is  permuted  into 
1  the  block  partitioned  form 

I 

r: 
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where  the  diagonal  blocks  are  square. 

In  our  case  we  use  the  following  dissection  strategy  in  order  to 
permute  the  observation  matrix  A  into  the  partitioned  form  (1.4) 

Our  procedure  will  be  called  nested  bisection. 

Given  a  geodetical  position  network  on  a  geographical  region  E  , 
first  pick  a  latitude  so  that  approximately  one-half  of  all  the  stations 
lie  south  of  this  latitude.  This  forms  two  blocks  of  interior  stations 
and  one  block  of  separator  or  junction  stations  and  contributes  one  level 
of  nested  bisection  (see  Figure  3). 


interior  stations 

junction  stations 
interior  stations 


FIGURE  3 


One  level  of  nested  bisection. 


Now  order  the  stations  in  ]R  so  that  those  in  the  interior  regions 
^  appear  first/  those  in  the  interior  region appear  second,  and 
those  in  the  junction  regionJB  appear  last;  order  the  observations 
(i.e.,  order  the  equations),  so  that  those  involving  stations  in 


come  first,  followed  by  those  involving  stations  in«/l»2  ;  then  the 

observation  matrix  A  can  be  assembled  into  the  block- partitioned  form: 

m  m 

A  = 


Thus  A  can  be  expressed  in  the  block-partitioned  form: 


where  the  A^  contains  nonzero  components  of  equations  corresponding 
to  coordinates  of  the  interior  stations  in«^^  and  where  the  contain 
the  nonzero  components  of  equations  corresponding  to  the  coordinates  of 
the  stations  in  the  junction  region  JB  . 

Next,  in  each  of  these  halves  we  pick  a  longitude  so  that  approximately 
one-half  of  the  stations  in  that  region  lie  to  the  east  of  that  longitude. 

This  constitutes  level  2  of  nested  bisection.  The  process  can  then  be 
continued  by  successively  subdividing  the  smaller  regions,  alternating  between 
latitudinal  and  longitudinal  dividing  lines.  Figure  k  illustrates  three  levels 
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It  follows  that  if  nested  bisection  is  carried  out  to  k  levels, 
then  the  partitioned  form  of  the  assembled  observation  matrix  has: 

i)  2  diagonal  blocks  associated  with 
interior  regions,  and 

ii)  2^-1  blocks  associated  with  junction  regions. 

In  particular,  there  are 

iii)  2  ”  junction  blocks  which  are  each  coupled  to 
2  interior  regions,  said 

\  k_i 

iv;  2  -1  junction  blocks  which  are  each  coupled  to 

k  interior  regions. 

Heuristically,  one  normally  would  like  to  perform  the  bisection 
process  so  that  the  sets  of  junction  stations  are  minimal  at  each  level, 
thus  maximizing  the  numbers  of  columns  in  the  diagonal  blocks.  The  process 
is  stopped  at  the  level  k  at  which  the  2  diagonal  blocks  are  suffi¬ 
ciently  dense  or  at  the  level  at  which  further  subdivisions  are  not 
feasible  or  are  not  necessary  for  the  particular  adjustment  problem. 

Our  proposed  block  orthogonal  decomposition  algorithm  for  an  obser¬ 
vation  matrix  A  already  in  the  partitioned  form  determined  by  nested 
bisection  is  deferred  to  the  next  section. 

3-  The  Block  Orthogonal  Decomposition. 

In  this  section  we  describe  a  block  orthogonal  decomposition  algorithm 
for  solving  the  least  squares  adjustment  problem  min||b-Ax||g  ,  where 
the  observation  matrix  A  has  been  assembled  into  the  general  block 
diagonal  form  (l.U).  Here  we  assume  that  the  structure  of  A  is  specified 
by  the  nested  bisection  scheme  described  in  Section  2.  Other  dissection 


schemes  may  he  preferable  in  certain  applications  (see  Golub  and 
Plemmons  [1980]). 

We  first  illustrate  the  method  with  k  =  2  levels  of  nested 
bisection,  as  given  in  Figure  5* 


FIGURE  5 

Two  levels  of  nested  bisection. 


Suppose  that  the  associated  observation  matrix  A  is  assembled  into  the 
corresponding  block- partitioned  form,  giving 

B1  Dl" 

B2  °2 

C3  D3 

\  C4  \ 

Then  by  the  use  of  orthogonalization  techniques  based  upon,  for  example. 
Householder  reflections,  Givens  rotations  or  modified  Gram-Schmidt  ortho 
gonalization,  the  reduction  of  A  to  upper  triangular  form  proceeds  as 


follows: 


At  the  first  stage,  each  diagonal  block  is  reduced  by 
orthogonal  transformations. 


Here  the  Q.  are  orthogonal  matrices  (which  of  course  need  not 

1  t  T-y 


t  f  Ri  ' 

be  formed  explicitly)  and  A^  =  , 


where  R^,  = 


,  yielding 


The  row  blocks  corresponding  to  the  upper  triangular  matrices 
are  then  merged  through  a  permutation  of  the  rows,  yielding 


This  completes  the  first  stage  of  the  reduction.  For  the  intermediate 


stages,  pairs  of  merged  blocks  corresponding  to  junction  stations  are 


reduced.  First , 


are  reduced  to  upper  triangular 


form  by  orthogonal  transformations,  yielding 
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Sit mm 


hi  i*— rinii 


Then  merging  the  triangular  factors  R^  and  Rg  through  a  permutation 
of  the  rows,  yields 


To  complete  the  intermediate  stages, 
form  by  orthogonal  transformations. 


is  reduced  to  upper  triangular 
yielding 


-4  °4 


v  > 


Here  R  is  the  Cholesky  factor  for  A  .  Let  n^  denote  the 
order  of  R^  for  1  =  1,  ...,7  and  let  =  (e^, . . . ,Cy)^  denote  the 
result  of  applying  the  same  sequence  of  orthogonal  transformations  and 
permutations  to  b  ,  where  each  c^  is  an  n^-vector.  For  the  final 
step  of  the  solution  process,  the  least  squares  solution  y  to  Ax  b 
is  computed  as  follows. 

Partition  y  as  y*  =  (y1,...,y^)t  where  yi  is  an  ^-vector, 
i  =  1,...,7.  Then  the  following  upper  triangular  systems  are  solved 
successively  by  back- substitution  for  the  vectors  ,  i  =  7,6,...,1  . 

Ryyy  -  Crj  f 

=  ci  "  >  i  =  6>5> 

Riyi  =  ci  “  Ciy6  "  Diy7  '  1  =  4,5, 

Riyi  =  ci  -  b^5  ”  D°y7  »  i  =  2,1. 

The  general  reduction  process  is  described  next  in  terms  of  three 
basic  steps.  Let  A  and  b  denote  the  observation  matrix  and  vector 
resulting  from  k  levels  of  nested  bisection  of  a  geodetic  position 
network  on  some  geographical  region.  Assume  that  A  has  been  assembled 

into  the  general  block- partitioned  form  (l.U),  with  2k  diagonal  blocks 

It  K 

and  2  -1  remaining  column  blocks.  Letting  t  =  2  ,  we  write  A  as 
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Al,t+1  *  '  *  Al,2t-1 

A2,t+1  '  '  '  A2,2t-1 


A  = 


(3.1) 


t,t+l 


t,2t-l 


For  a  certain  flexibility  of  the  algorithm  and  also  for  simplicity  in  the 
notation,  we  do  not  altogether  distinguish  here  between  zero  and  nonzero 
blocks  A.  .  .  The  zero  pattern  of  these  blocks  depends  on  the  number  of 

levels,  k  ,  to  which  the  nested  bisection  process  is  carried.  Particular 
attention  was  paid  to  this  pattern  for  the  case  of  k  =  2  levels  in  the 
illustrative  example  just  completed. 


Algorithm  1.  This  algorithm  computes  the  Cholesky  factor  R  and  the  least 
squares  solution  y  to  Ax  w  b  where  A  results  from  k  levels  of 
nested  bisection  and  A  has  the  block  form  (3-l),  with  t  =  2k  . 


Step  1.  Reduce  each  diagonal  block  A^  of  A  to  upper  triangular  form 
by  orthogonal  transformations  and  merge  the  reduced  blocks, 
l)  Do  for  i  =  1,2,..., t. 


l)  Determine 
(Note  that 


so  that 


❖i 


0  * 

QT  need  not  be  formed  explicitly). 
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2 )  Compute 


WAi,t+l'“*'Ai,2t-l] 


R  A0  a0 

I  Ai,t+l"*  Ai,2t-1 


0  A.  ,  -  •  •  *  A.  ,  « 
i,t+l  i,2t-l 


2)  Merge  the  reduced  row  blocks  by  row  permutations  so  that  the 
resulting  matrix  has  the  form 


1, t+l  l,t+2 

\°  A0 

2, t+l  2,t+2 


t,t+l 

S.,t+1 

^2,t+l 


t ,  t+2 


3,t+2 

%,t+2 


Al,3t/2  Al,(3t/2)+l 

A2,3t/2  A2,(3t/2)+l 


A°  / 
t,3t/2 


At-l,3t/2 

At,3t/2 


At,(3t/2)+l*  ' 

Al, (3t/2)+l  *  ' 

A2,(3t/2)+l  '  * 

A3,(3t/2)+l  •  * 

\,(3t/2)+l  '  * 


S;-l,  (3t/2)+l 

^t,(3t/2)+l  * 


l,2t-l 

1° 

2,2t-l 


0 

t,2t-l 

1 

l,2t-l 

1 

2,2t-l 

1 

3»2t-l 

V,2t-1 


1 

t-l,2t-^ 
1 

t,2t-l 
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Step  2.  Reduce  and  merge  the  intermediate- stage  blocks, 
l)  Do  for  u  =  t,  t/2,...,  t/2k_1  =  2  . 

l)  Do  for  v  =  l,3,...,u-l 

I  l)  Reduce  each  pair  of  row  diagonal  blocks 


to  upper  triangular  form  by  orthogonal  transformation, 
as  in  Step  1. 

2)  Merge  the  resulting  reduced  row  blocks  by  row  permutations 
so  that  the  upper  triangular  blocks  R^  appear  first,  as 
in  Step  1. 

At  the  end  of  Step  2,  A  has  been  reduced  by  orthogonal  transformations 
to  the  following  form,  where  each  R^  is  upper  triangular  and  where  certain 
of  the  blocks  A?,  are  zero. 

ij 
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The  reduction  algorithm  just  described  for  the  observation  matrix  A 

can  be  interpreted  from  a  network-reduction  viewpoint  as  follows.  Suppose 

that  A  results  from  a  nested  bisection  of  the  geographical  region  to  k 

levels.  Then  at  the  first  stage  of  the  reduction  process,  orthogonal 

k 

transformations  are  applied  to  each  of  the  2  blocks  corresponding  to 
the  interior  regions ,  to  reduce  the  coordinates  of  stations  not  coupled 
to  stations  outside  that  block  by  an  observation.  Modified  junction  stations 
in  the  separator  blocks  are  kept  until  nearby  interior  blocks  are  reduced. 
Then  clusters  of  blocks  of  junction  stations  are  grouped  together  (merged) 
to  form  higher  level  blocks.  At  the  intermediate  stages  of  the  reduction 
process,  some  station  coordinates  are  now  interior  and  can  be  reduced  by 
orthogonal  transformations.  The  process  continues  until  at  the  last  stage 
the  remaining  stations  are  all  interior  and  their  coordinates  can  be  reduced. 
At  this  point  A  is  completely  reduced  by  orthogonal  transformations  to  its 
Cholesky  factor  R  ,  and  correspondingly,  the  vector  b  is  reduced  to  c 
as  indicated  in  Step  3.  To  determine  the  least  squares  solution  y  to 
Ax  sa  b  ,  the  process  is,  in  a  sense,  reversed  to  back  substitute  the  co¬ 
ordinates  to  successively  lower  levels  until  all  of  the  corrections  have 
been  found. 

Notice  that  at  each  stage  of  the  reduction  process  it  is  possible  to 
obtain  a  "diagnostic  solution"  (see  Meissl  [1979]).  Here  we  hold  the  co¬ 
ordinates  of  the  junction  stations  fixed  and  solve  for  the  coordinates  of 
the  reduced  interior  stations  at  that  stage. 

We  emphasize  again  that,  fora  certain  flexibility,  full  advantage  has 


not  been  taken  in  Algorithm  1  of  the  zero  pattern  of  the  blocks 

of  A  as  given  by  (3.l).  This  pattern  of  course  determines  the  block 

structure  of  the  Cholesky  factor  R  of  A  as  given  by  (3.2).  Basically, 

k+1 

R  has  the  same  type  of  block  structure  as  A  ,  but  with  2  -1  upper  - 

triangular  diagonal  blocks.  For  nested  bisection  to  k  =  k  levels, 
where  A  is  assembled  into  the  form  (2.3),  the  Cholesky  factor  R  has 
the  following  structure . 


In  order  to  facilitate  an  analysis  of  the  results  of  a  least 
squares  adjustment,  it  is  often  desirable  to  compute  some  or  all  of 
the  elements  of  the  variance-covariance  matrix  (A^A) .  Since 


oAo-1  =  (r^)-1  , 


the  special  block  structure  of  R  just  discussed  can  be  used  advanta¬ 
geously  in  computing  the  variances  and  covariances.  Such  a  procedure 
is  given  in  the  next  section  for  a  more  generally  sparse  Cholesky  factor  R  . 

4.  Computation  of  the  Variances. 

In  many  adjustment  problems  (see,  for  example,  Hanson  [1978])  it  is 
necessary  to  compute  the  variances  and  covariances  associated  with  the 
regression  coefficients  in  order  to  estimate  the  accuracy  of  the  results. 
Under  the  usual  assumptions,  the  variance  of  the  i-th  coefficient  is  pro¬ 
portional  to  the  (i,i)  element  of  (AtA)"1  .  If  R  is  sparse,  then  the 
diagonal  elements  of  (A^A)”1,  can  be  calculated  quite  efficiently.  Indeed, 
it  is  easy  to  compute  all  the  elements  of  (A^A)-1  which  are  associated 
with  the  non-zero  elements  of  R  ,  the  Cholesky  factor.  We  describe  the 
procedure  next. 

Using  the  orthogonalization  algorithm  we  determine  the  Cholesky 
factor  R  so  that 


AtA  =  RtR  . 


Suppose 


when  (i, j )e  K 
when  (i,jH  K  . 
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ill  fi  iF 


Our  objective  is  to  determine 

((AtA)“1}ij  When  (i,j)  e  K  . 

Let  us  write 

(A^)-1  =  Z  =  [z1>...,zn]  , 

where  z^  is  the  i-th  column  of  the  matrix  Z  . 

Since 

A^Z  =  I 
RZ  =  (Rt)“1  . 

Note  that 

■ 1  /  -n  • 

From  (4.l)  and  (4.2),  we  see  that 

R  zn  *  \  *  (rnnrl  <4  ‘  0,1> 

so  that  we  can  solve  for  by  back  substitution.  Thus 

z  =  (r  )"^ 

nn  nn 

and  for  i  =  n-l,n-2,. . . ,1 


(4.1) 


(4.2) 


in 


-I  i 

j=i+l 


ii 


in 


-I 


J-i+1 

(i,j)e  K 


lii 

rii 


in 
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t-’ 


•  ? 

I 


i  I 


r: 


Let  I  =  min  fi  |  r.  ^  0}  .  It  is  possible  to  calculate 
1  <  i  <  n-1  in 

z^n  for  i  =  n-l,n-2, . . . ,In  .  Once  these  components  have  been  computed, 
it  is  only  necessary  to  save  those  elements  for  which  (i,n)  €  K  . 

Note 


z  =  z 
in  m 


Now  assume  we  have  calculated  those  elements  of  z  ,z  ...... z.  , 

n7  n-17  7  f +1 


for  which 


r  ^  °  When  p=J-7'"'n  >  q  =  *+•!,..,  n  . 


Thus,  by  symmetry  we  have  computed 


z  .  for  q  >  l  and  (f,q)e  K  . 
ql 


Now  for  i  =  1,2,.,  1-1 


and 


I 

j=i 

n 


r .  .  z .  =  0 

ij  01 


V*  1 

L.  |.j  3t  T. . 


S=t 


Hence 


n 


=  -i-  (  jl  _y 

Tn  Tu  *— 


n 


r  .  z.  ) 
LA  31 


A=*+ 1 


=  -i-  (— i-  -  V*  r  .  z.  ) 
ra  r«  L  lS  d' 


J-i+ 1 
(i,j)€  K 
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• 


T 


Let  I.  =  min  {i  |  r.  ,  £  0}  .  Then  for  i  =  4-1,.,., I 
1  1  <  i  <  4-1  1  1 


“if 


d-i+i 

(i,j  )  €  K 


ij  jf 
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-I 


j=4+l 
(i,j)e  K 


/ru 


Again,  after  this  calculation  is  performed,  we  save  only  those  elements 
for  which  (i,4)  e  K  .  The  above  algorithm  thus  describes  a  method  for 
computing  the  elements  of  the  inverse  of  (A^A)  which  are  associated 
with  the  non-zero  elements  of  R  .  Such  a  procedure  can  be  quite 
efficient  when  compared  to  computing 

(A^)-1  =  R_1(Rt)'1  . 

For  example,  suppose  we  need  the  diagonal  elements  of  (A^A)”1  when 

r^  i  0  for  i  =  J  and  j  =  i+1  ,  and 
r .  .  =  o  otherwise, 

i.e.  R  is  bi-diagonal.  The  matrix  R_1  will  be  completely  filled  in 
above  the  diagonal  and  hence  0(n  )  numerical  operations  are  required  to 
compute  the  diagonal  elements  of  (A^A)'1  .  The  algorithm  we  have  outlined 

above  would  require  0(n)  operations.  Even  greater  savings  can  be  expected 
for  the  Cholesky  factor  R  of  the  form  (3.2),  resulting  from  nested  bi¬ 
section. 


5.  Final  Remarks. 

To  summarize,  an  alternative  has  been  provided  here  to  the  formation 
and  solution  of  the  normal  equations  in  least  squares  adjustment  problems. 

In  particular,  it  has  been  shown  how  a  block- orthogonal  decomposition  method 
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can  be  used  in  conjunction  with  a  nested  dissection  scheme  to  provide 
a  least  squares  algorithm  for  certain  geodetic  adjustment  problems. 

Some  well-known  advantages  of  dissection  schemes  for  sparse  linear  systems 
are  that  they  facilitate  efficient  data  management  techniques,  they  allow 
for  the  use  of  packaged  matrix  decomposition  routines  for  the  dense 
component  parts  of  the  problem, and  they  can  allow  for  the  use  of  parallel 
processing.  In  the  past,  the  combination  of  the  normal  equations  approach 
with  these  dissection  techniques  (in  particular  Helmert  blocking)  has  been 
preferred,  partly  because  of  tradition  and  partly  because  of  the  simplicity 
and  numerical  efficiency  of  the  Cholesky  decomposition  method.  However, 
the  use  of  an  orthogonal  decomposition  scheme  applied  directly  to  an 
observation  matrix  A  which  has  also  been  partitioned  by  a  dissection 
scheme  has  several  advantages  over  the  normal  equations  approach.  First, 
the  QR  orthogonal  decomposition  of  A  allows  for  an  efficient  and 
stable  method  o"  adding  observations  to  the  data  (See  Gill,  Golub,  Murray 
and  Saunders  [1974]).  Such  methods  are  crucial  in  certain  large-scale 
adjustment  problems  (see  Hanson  [1978]).  Secondly,  possible  fill-in  that 
can  occur  in  forming  the  normal  equation  matrix  A^A  is  avoided.  A 
statistical  study  of  such  fill-in  is  provided  by  Bjorck  [1976].  Meissl 
[1979]  reports  that  sane  fill-in  can  be  expected  in  forming  A^A  in  the 
readjustment  of  the  North  American  Datum.  This  problem  cannot  be  over¬ 
emphasized  in  such  large  scale-  systems  (6,000,000  equations  and  400,000 
unknowns).  But  perhaps  the  most  crucial  advantage  of  the  use  of  ortho¬ 
gonal-decomposition  schemes  here  is  that  they  may  reduce  the  effects  of 
ill-conditioning  in  adjustment  calculations. 
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In  this  paper  we  have  treated  only  one  aspect  of  nested  dissection 
in  least  squares  problems,  that  of  decomposing  a  geodetical  position 
network  by  the  process  of  nested  bisection.  However,  the  block  diagonal 
form  of  the  matrix  in  (l.4)  can  arise  in  other  dissection  schemes  such 
as  one-way  dissection  (see  George,  Poole  and  Voight  [1978]  for  a  description 
of  this  scheme  for  solving  the  normal  equations  associated  with  finite 
element  problems).  The  form  also  arises  in  other  contexts,  such  as  photo- 
grammetry  (See  Golub,  Luk  and  Pagano  [1979]).  Least  squares  schemes  based 
in  part  upon  block  iterative  methods  (see  Plemmons  [1979 Dor  a  combination 
of  direct  and  iterative  methods  may  be  preferable  in  some  applications. 
Moreover,  the  general  problem  of  permuting  A  into  the  form  (1.4 )  by 
some  graph-theoretic  algorithm  for  ordering  the  rows  and  columns  of  A 
(see  Weil  and  Kettler  [1971])  has  not  been  considered  in  this  paper. 

Some  of  these  topics  will  be  addressed  further  in  Golub  and  Plemmons  [1980]. 
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